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ABSTRACT 

We carry out a numerical hydrodynamical modeling for the evolution of a 
relativistic collimated outflow, as it interacts with the surrounding medium, and 
calculate the light-curve resulting from synchrotron emission of the shocked fluid. 
The hydrodynamic equations are reduced to 1-D by assuming axial symmetry 
and integrating over the radial profile of the flow, thus considerably reducing 
the computation time. We present results for a number of different initial jet 
structures, including several different power-laws and a Gaussian profile for the 
dependence of the energy per unit solid angle, e, and the Lorentz factor, T, on 
the angle from the jet symmetry axis. Our choice of parameters for the various 
calculations is motivated by the current knowledge of relativistic outflows from 
gamma-ray bursts and the observed afterglow light-curves. Comparison of the 
light curves for different jet profiles with GRB afterglow observations provides 
constraints on the jet structure. One of the main results we find is that the 
transverse fluid velocity in the comoving frame (vt) and the speed of sideways 
expansion, for smooth jet profiles, is typically much smaller than the speed of 
sound (c s ) throughout much of the evolution of the jet; v t approaches c s when T 
along the jet axis becomes of order a few (for large angular gradient of e, v t ~ c s 
while T is still large). This result suggests that the dynamics of relativistic 
structured jets may be reasonably described by a simple analytic model where e 
is independent of time, as long as Y along the jet-axis is larger than a few. 



Subject headings: gamma-rays: bursts - gamma-rays: theory 
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1. Introduction 

The great advance in our understanding of Gamma-Ray Bursts (GRBs) in the last 
five years has largely resulted from the observation and modeling of afterglow radiation - 
emission observed for days to months after the end of a GRB, in the X-ray, optical and 
radio bands. The basic procedure for obtaining information about the explosion, such as 
the energy release, opening angle of the emergent jet, the density of the medium in the 
immediate vicinity of the GRB etc., is by comparing the observed afterglow light-curve 
with the theoretically calculated flux (Wijers & Galama 1999; Granot, Piran & Sari 1999b; 
Chevalier & Li 2000; Panaitescu & Kumar 2001a,b, 2002). Most works on GRB jets assume 
a homogeneous (or 'top hat') jet, where all the hydrodynamic quantities of the jet, such as 
its Lorentz factor and energy density, are the same within some finite, well defined, opening 
angle around the jet axis, and drop to zero at larger angles. 

A comparison of theoretically calculated light-curves, under several simplifying assump- 
tions described below (and assuming a 'top hat' jet), with observed light-curves in X-ray, 
optical, and radio bands for 8 GRBs, has led to a number of interesting results (Panaitescu 
& Kumar 2001b). Perhaps the most remarkable discovery is that the kinetic energy in the 
relativistic outflow is nearly the same, within a factor of 5, for the set of eight GRBs. A 
similar result has been obtained by Piran et al. (2001) through a method which requires 
fewer assumptions. Frail et al. (2001) have also found that the energy radiated in GRBs 
does not vary much from one burst to another. The opening angle for GRB jets is found to 
be in the range of 2-20 degrees, and the density of the external medium in the vicinity of 
GRBs is estimated to be between 10~ 3 and 30 cm -3 . Moreover, there is no firm evidence 
for the density to vary as inverse squared distance in all but one case (Price et al. 2002; 
Panaitescu & Kumar, 2002), which is surprising in light of the currently popular model for 
GRBs - the collapsar model. 

The possibility that GRB jets can display an angular structure, i.e. that the Lorentz 
factor, 7, and energy per unit solid angle, e, in the GRB outflow can vary smoothly as power 
laws in the angle 6 from the jet axis, was proposed by Meszaros, Rees & Wijers (1998). 
Recently, in view of the evidence described above for a roughly constant energy in the 
gamma-ray emission and in the kinetic energy of the afterglow shock, it has been suggested 
that GRB jets might have a universal structure, and the differences in the observed properties 
of GRBs and their afterglows arise due to different viewing angles, 9 ohs , w.r.t the jet axis 
(Lipunov, Postnov & Prokhorov 2001; Rossi, Lazzati & Rees 2002; Zhang & Meszaros 2002). 
In this interpretation, the jet break in the light curve occurs when the Lorentz factor along 
the line of sight, 7(# bs), drops to ~ 9~^ s , so that the jet break time, tj, is determined by 
the viewing angle, # b s , rather than by the opening angle of the 'top hat' jet, as in the 



- 3- 



conventional interpretation. 

The calculation of light-curves from a shock-heated, collimated, relativistic outflow has 
been carried out by a few research groups (Rhoads 1999; Panaitescu & Meszaros 1999; Kumar 
& Panaitescu 2000; Moderski, Sikora & Bulik 2000; Granot et al. 2002). However, most 
of the works to date have been based on a simplified model for the jet dynamics and on a 
number of ad-hoc assumptions. All the above works assume a 'top hat' jet, and furthermore, 
most of them model the dynamics of the jet at times much greater than the deceleration time 
as uniform expansion at the sound speed or the speed of light (in the local rest frame of the 
shocked fluid)- the results are nearly the same for both of these cases. Similar simplifications 
were made in the recent work on a universal structured jet (Rossi, Lazzati & Rees 2002). An 
exception to this, is the work of Granot et al. (2001), where the dynamics of an initial 'top 
hat' jet were calculated using a hydrodynamic simulation, and the resulting light curves were 
calculated numerically. However, such hydrodynamic simulations are very time consuming, 
and difficult to apply to a structured jet, so that there is currently no rigorous treatment of 
the hydrodynamic evolution of a structured jet. In this paper we develop such a rigorous 
treatment for the dynamics of structured jets, which at the same time is not very time 
consuming and may become practical to include in fits to afterglow observations. 

Another simplification made in previous works (including all the works mentioned 
above), and in the lack of a better alternative, is also made in this work, is that the strength 
of the magnetic field and the energy in the electrons are determined by assuming that the en- 
ergy densities of the magnetic field and of the electrons are constant fractions of the internal 
energy density of the shocked fluid. It is unclear how some of the simplifying assumptions 
in the afterglow light-curve modeling effect the overall burst parameters and properties we 
have inferred as described above. 

Some progress has been made recently toward understanding the generation of magnetic 
fields in relativistic collisionless shocks: the numerical simulations of Medvedev (2002) show 
that magnetic fields generated behind collisionless relativistic shocks via the Weibel instabil- 
ity (Medvedev & Loeb 1999) do not decay to very low values within a short distance behind 
the shock, as was previously thought (Gruzinov 1999, 2001), but rather approach a finite 
value in the bulk of the shocked fluid behind the shock, which might be compatible with the 
values inferred from afterglow observations. Moreover, the modeling of GRB afterglow light- 
curves indicates that the energy fraction in electrons is close to equipartition (Panaitescu 
& Kumar, 2001b), hence the parametrization of electron energy does not appear to be a 
serious drawback for current models. Thus, at present, one of the biggest uncertainties in 
the afterglow modeling is the assumption of a uniform jet and the simplified jet dynamics. 
The purpose of this paper is to remedy this situation and develop a much more realistic 
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model for GRB jets. Fitting afterglow observations with light-curves that are calculated 
using a realistic jet model & dynamics may both constrain the structure of GRB jets (the 
initial distribution of the Lorentz factor and energy per unit solid angle as a function of the 
angle from the jet axis), and provide more accurate estimates for the physical parameters, 
which include the external density profile and the parameters describing the micro-physics 
of relativistic collissionless shocks. 

In the next section (§2) we discuss the evolution of structured jets. In §2.1 we describe 
our hydrodynamical scheme where we begin from the full hydrodynamic equations, assume 
axial symmetry and integrate over the radial structure, thus reducing the problem to a set 
of one dimensional partial differential equations that are solved numerically. The initial and 
boundary conditions are outlined in §2.2, while results for some physically interesting cases 
are shown in §2.3. In §3 we describe the light-curve calculation and compare the results of 
hydro simulations with a simplified model. The main results are summarized in §4. 



2. Jet modeling 

We begin with a brief description of the uniform jet model, and then we describe in 
some detail the evolution of a more realistic, structured, jet and the afterglow light-curves 
resulting from emission by the shock heated medium swept up by the jet. 

Most calculations of GRB light-curves have assumed that the properties of the relativis- 
tic outflow do not vary across the jet, and that the jet dynamics is described by a uniform 
lateral expansion in the comoving frame, at close to the speed of sound, c s , which for a hot 
relativistic plasma is 3 -1 / 2 times the speed of light, c. These assumptions drastically simplify 
the calculation of the evolution of the jet opening angle, 9j, with time: the increase in the 
lateral size of the jet in comoving time 5t co is c s 5t co , and so the change to its angular size is 
50 j = c s 5t co /r = (c s /c)5r/(r^y), or, 

d9j c s 
dr cryr 

This equation, together with the energy conservation equation, describe the dynamics of a 
uniform relativistic disk or a jet. The implication of this equation is that the jet opening angle 
9j starts to increase when 7 drops below 0J 1 , and from that time onward the jet opening angle 
is roughly 7 -1 . A detailed discussion of the uniform jet dynamics and lightcurve calculation 
can be found in a number of papers (e.g. Rhoads 1999; Panaitescu & Meszaros 1999; Sari 
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et al. 1999; Kumar & Panaitescu 2000). Such a uniform jet with sharp well defined edges 
shall be referred to as a 'top hat' jet. 

However, 2D hydrodynamical simulations of the evolution of a jet that is initially uni- 
form within some finite opening angle (Granot et al. 2001) have shown that the lateral 
expansion of the jet is smaller than that predicted by the simple models described above. 
This suggests that the assumption of lateral expansion at close to the sound speed in the 
comoving frame, that is made in most simple jet models, is not valid. Nevertheless, the light 
curves calculated from these simulations show a sharp jet break in the light curves, similar 
to that seen in most afterglow observations, around the time 7 drops to Qj 1 . 



2.1. Dynamics of Structured Relativistic jets 

Clearly, it is unrealistic to assume that the outflow from GRB explosions will be uniform 
within some finite opening angle, outside of which the Lorentz factor, 7, and energy per unit 
solid angle, e, decrease very sharply (i.e. a 'top hat' jet). A more realistic situation is that the 
Lorentz factor (LF), the energy density etc. are smooth functions of the angle, 8, from the jet 
axis, and possibly also of the distance, r, from the central source. The causality consideration 
suggests that the outflow is unlikely to be uniform over large angles, and moreover it provides 
a limit on how rapidly initial inhomogeneities can be smoothed out. Let the LF of the shell 
after elapsed time t since the explosion, measured in the lab frame, be 7(t). The comoving 
time corresponding to this is ~ t/7, and the distance traversed by sound waves during this 
interval is c s t/^f ~ ct/3 1 ^ 2 7. Therefore, the angular size of a causally connected region is 
~ l/3 1 / 2 7, and inhomogeneities on an angular scale of 9^ > 7 -1 , if present initially, will 
persist; the inhomogeneities can be smoothed out only when the LF has dropped below 9~^. 
As an example, the large angular scale inhomogeneities for a jet of opening angle 5° start 
to decrease only when the bulk LF has dropped below ~ 10, or roughly one day after the 
explosion (as seen by the observer). It should also be noted that if one were to start with a 
uniform jet, or a top-hat profile for the LF or e, the large gradient at the edge will decrease 
with time and the jet will develop angular structure (e.g. Granot et al. 2001). 

The remainder of this section is devoted to the solution of the relativistic hydrodynamic 
equations to describe the evolution of jets from GRBs. The starting point is the relativistic 
fluid equations (e.g. Landau & Lifshitz, 1959): 

7% = , T^ = wuV + pg^ , (1) 

where is the energy-momentum tensor for an ideal fluid, is the 4- velocity of the fluid, 
g^ v is the metric tensor, p is the pressure and w = p + e + p is the proper enthalpy density, 
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where p and e are the proper rest mass density and internal energy density respectively, 
and c = 1 in our units. We use a spherical coordinate system and assume the flow possesses 
axial symmetry about the z-axis, i.e. u^,d/d(j) = 0. Under these assumptions the t,r and 9 
components of equation (1) are 

^-(wT 2 -p) + l|-(r 2 W rV) + —^-:-^(sm9wT 2 v e ) = , (2) 
ot r 2 or r sin 9 09 

9 I r2 \ , 1 9 i 2 F 2 i\ . 1 d . 2 dp wY 2 v 2 e 

— (wT v r ) + -y-~-(r wT v r ) H — — {sm9wT v r v e ) + -= = , (3) 

at r 2 or r sin 9 09 or r 

9, rf , 1 d , 2 ^ 2 . 1 9 , . . rf I dp wT 2 v r v e 

T-(wT 2 v e ) + - — (r 2 wT 2 v r v 6 ) + —- — ( sm 9wT 2 v 2 ) + + = , 4 

ot r 2 Or r sin 9 09 r 09 r 

where v r and v e are the r and 9 components of the fluid velocity, and r = (1 — v 2 — Vq)^ 1 ^ 2 
is the Lorentz factor of the fluid. Assuming that pair production has a negligible effect on 
the rest mass density, baryon number conservation implies 

We assume an equation of state 

7-1 

P = (7 - l)e = — — - p) 

7 

Equations (2)-(6) can be solved together to determine the structure and evolution of 
the outflow from GRBs. The computation time, for a 1 GHz clock speed computer, and 
for a modest resolution in r & 9 coordinate of 100x1000 (in order to keep the error small 
in finite difference schemes) and 5000 time steps, is expected to take of order several hours 
to complete one run for one set of initial conditions; for comparison the 2-D relativistic jet 
hydrodynamics calculation of Miller and Hughes, reported in Granot et al. (2001), took 
several hours to days of computation time, for low to medium resolution runs, to follow the 
evolution for > 10 observer days, while an even longer computational time was required 
for the higher resolution runs. The successful modeling of light-curves of a single GRB 
to determine various parameters requires several thousand runs, and thus the computation 
time to model one GRB, using a 2-D code, is currently estimated to range between months 
to years. Using many processors simultaneously can help reduce the actual overall time 
required, but at any rate, this requires a great computational effort. 



with 



7 



4r + 1 
3r 



(6) 
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The computation time can be drastically decreased by reducing the problem to a 1-D 
system, by integrating out the radial dependence for all of the relevant variables, over the 
width of the outflow plus the swept up material. The physical motivation for this is that 
jets in GRBs are in fact thin shells. 1 This procedure reduces the computing time drastically 
without introducing a significant loss of information as far as the emergent synchrotron 
emission is concerned; we find that the lightcurves from a relativistic spherical shock which 
has radial structure described by the self-similar Blandford-McKee (1976) solution is almost 
the same as in a model where the radial dependences have been integrated out and the shell 
thickness is taken to be zero (see Figure 5 of Granot, Piran & Sari 1999a). 

The shock front is a two dimensional surface described by r = R(6,t). The shocked 
fluid is concentrated in a thin shell of thickness AR ~ R/AT 2 <C R for a relativistic flow 2 , 
and therefore it makes sense to integrate all the dependent variables, such as p, w, T, etc., 
over the width of the shell in the radial direction. We define quantities averaged over r, at 
a fixed 9 and lab frame time t, as follows: 



pR _ pR _ pR 

n= / drr 2 p , nr 2 = / drr 2 pT 2 , Tu e U = / drr 2 pYu e , 
Jo Jo Jo 



pR pR pR 

nrw r = / drr 2 pTu r , Xi^-^rUe = / drr 2 pu r u e , X2ujjll = / drr 2 pu 2 B , (8) 
Jo Jo Jo 



(7) 



_dM s [« _dM f K 2 

^ = ^n=J drr pv ' ^ = ^ = 1 drrpoT > (9) 

where /i s and /x are the rest mass per unit solid angle of swept-up material and of the 
original ejecta, respectively, and xi & Xi are dimensionless correlation coefficients, of order 
unity magnitude, which are taken to be independent of time. Integration of equation (2) 
times r 2 , over the radial interval corresponding to the width of the shell, when the shell is 
located at R(t), yields 



d_ 
di 



er 2 - n 



1 d 



sin 9 ueTt; 
R 



PcAR)R 2 ^r , (10) 



*At a distance r from the center of the explosion, the laboratory frame radial thickness of the ejecta plus 
the swept-up shock heated material moving with LF 7 is ~ r/Aj 2 , whereas its transverse dimension is rOj. 
Therefore the geometric shape of the system is that of a thin disk as long as 0j ^> I/47 2 . 

2 In fact, even as the flow becomes non-relativistic, we still expect AR/R < 0.1, as in the Sedov- Taylor 
self-similar solution, so that the thin shell approximation should still be reasonable. 
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where 



r ' 7-1 

An integration of equations (3) and (4) over the width of the shell gives 



(11) 



9 r^-l 1 9 

at L J sm 9 89 



sin 9 Xi u eu r -i 
R 



tt2 



2IT 

n ~r =0. 



(12) 



^ r - n Id 
[true] + 



dt 



sin 9 89 



sin9x2 U 2 e C 
R 



+ i? + R 89 ~ ' 



(13) 



and the closure relation, given below, is obtained by integrating the equation T 2 = I+u^+Uq, 
times the pressure p, over the width of the shell 



r - Xsu 2 r - xiu\ = 1 • (14) 
Where %3 ~ 1 is a constant factor; X3 — 1 corresponds to the assumption that 1 — u r /T <^ 1. 



In deriving these equations it was assumed that the ejecta moves with the shocked ISM, 
i.e. that the radial and the transverse components of the ejecta velocity are same as those 
of the swept up ISM. Under this assumption the mass continuity equation for the ejecta and 
the shocked ISM are respectively 



dfio 1 



8 



dt sin 9 89 



sin 9 ue/io 
RT 



0, 



(15) 



dji s 1 d 



dt sin 9 89 



sin 9 uqHs 
RT 



. . d2 8R 
p cxt {R)R , 



(16) 



The velocity of the shock front is given by the shock jump conditions (Blandford & 
McKee 1976): 

= j , (2 . fl( r,.i)Tr 

V (r P , + 1) Wr,,, - 1) + 1] 2 

where T ps ps f is the post shock LF. The shock jump conditions imply that in the rest frame 
of the fluid before the shock (which in our case is the lab frame), the direction of the velocity 
of the fluid just behind the shock is always perpendicular to the shock front. In order to 
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propagate the shock front in time, we assume that the average velocity of the shocked fluid 
is a good approximation for its value just behind the shock, and obtain 

dR v sh 



(18) 



Equations (10)-(18) are solved numerically with appropriate initial conditions (discussed 
below), to determine the evolution of the jet. 



2.2. Initial & boundary conditions 

The initial conditions are chosen at a lab frame time to, after the internal shocks have 
ended, and before there is a significant deceleration due to the sweeping up of the external 
medium. We implement a number of different initial conditions, one of which is that the 
initial energy (including the rest mass energy) per unit solid angle, e, and the initial Lorentz 
factor (minus 1) are power-law functions of 9, outside of a core of opening angle C . The 
initial energy per unit solid angle, e(9,to), and the LF for the power-law model are 

e(0,f o ) = e e- a , r(0,to) = i + (r -i)e- fe , (19) 

where eo and r are the initial energy per unit solid angle and Lorentz factor at the jet axis, 
and 




Another initial condition we explore is a Gaussian profile for which e(6,to) and F(0,to) 
are proportional to exp(— 9 2 /29 2 c ). The Gaussian jet profile was mentioned in Zhang & 
Meszaros (2001), however they did not calculate the jet dynamics or lightcurve for this case. 

Equation (19), or its counterpart for the Gaussian case, are applied only as long as 
F(9,t ) > 1.1. At larger angles, we assume a uniform outflow, the parameters of which are 
set by the continuity condition. 

We assume that the velocity is initially purely in the radial direction, i.e., 

v e (9,t ) = , v r (9,t ) = ^l-T- 2 (9,t ) , (21) 

while the initial radius is given by R(9,to) = t v r (9,t ). 

The angular derivative of all dependent variables except ug is zero at the pole and the 
equator, whereas dug/ 89 at the pole is determined by the assumption of axisymmetry, and 
at the equator by the reflection symmetry; u e vanishes at the pole and at the equator. 
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2.3. Numerical Results 

Equations (10)-(18) are solved using the two step Lax-Wendroff scheme, to determine 
the evolution of the jet, for several different choices of initial conditions. The number of grid 
points in the angular directions is taken to be about 2000, and the time step is chosen to 
satisfy the Courant condition. The numerical solution respects global energy conservation 
to within 0.1%. We have looked into the dependence of the solution on the value of the 
dimensionless correlation parameters xi an d X2, and find consistent solutions for \i between 
0.8 and 1, and X2 approximately between 0.9 and 1.2. Outside of this range of parameters 
the code is unstable and the solution unphysical, and the energy conservation is not satisfied. 
Inside this range the solution is not sensitive to the exact value of Xi an d X2- 

The evolution of r(0), u$(8), and e(6) are shown in figures 1-4 for (a, b) = (0,2), (2,0), 
(2,2), and for a Gaussian jet. Note that the transverse velocity in the comoving frame of 
the shocked fluid, v' e = Ug, is much less than the sound speed, 3 -1 / 2 , throughout much of 
the time, and approaches the sound speed only when the jet Lorentz factor on the axis has 
fallen to a value of order a few. Clearly, this result depends on the gradient of the initial LF 
or e at the initial time, and therefore the transverse velocity is found to be largest for the 
Gaussian case, which has the highest gradient of all the models we have considered. 

The small value for v' g can be understood from equation (13). Ignoring the second order 
term in ug, and noticing that the "source term" for ug is the gradient of II, we find that 
ug ~ (T56)" 1 , where 56 is the angular scale for the variation of II or the energy density. 
Thus, we get an appreciable transverse velocity in the comoving frame only when f 56 < 1. 

Another way of deriving this result is from the shock jump conditions, that imply that 
the velocity, v ps , of the fluid just behind the shock in the rest frame of the fluid before the 
shock (the lab frame in our case) is perpendicular to the shock front. This implies that the 
angle a between v ps and f (i.e. v ps • f — cos a) satisfies 

v e 1 dR . . 

t!ma = vr-Rm- (22) 

For a relativistic flow vg < <C 1 so that v g <C v r ~ 1 and a ~ tana pa vg. Thus we have 
v e « -dlnR/dd, and since R « (1- l/2r 2 )t, this implies v e « -T^dY/dd, or v e ~ l/(rW) 
and ug = Tvg ~ 1/(T58), where 56 is the angle over which T varies appreciably. Therefore, 
this reproduces the result of the previous paragraph, as it is easy to show that the angular 
scale for the variation of II and T are similar. From the definition of II (equation 7), we have 
II ~ pR 2 AR ~ pR 3 /AT 2 and since the shock jump conditions imply 3p = e = AT 2 p ext (R), 
this gives II ~ p ext (R)R 3 /3, i.e. II is a power law in R w (1 — l/2r 2 )t, and therefore varies 
over the same angular scales as T. 
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The transfer of energy from small to large angles over the course of the evolution of the 
jet, from highly relativistic to mildly relativistic regime, is also found to be small (Fig. 3). 

The Gaussian initial jet profile is a smooth and more realistic version of a 'top hat' 
jet, where the hydrodynamic quantities are roughly constant within some typical opening 
angle, and sharply drop outside of this angle (though they drop smoothly, and not in a step 
function as in the 'top hat' jet). We therefore expect the results for an initial Gaussian 
profile to be similar to those of an initial 'top hat' jet profile. The hydrodynamic evolution 
of the latter has been investigated using a 2D hydrodynamic simulation (Granot et al. 2001) 
and was found to be quite similar to our results, namely the lateral spreading of the jet was 
much smaller than the prediction of simple 'top hat' jet models, and there was very little 
lateral transfer of energy. The fact that our hydrodynamic results for an initial Gaussian 
profile are similar to the results of 2D hydrodynamic simulations of an initial 'top hat' jet 
profile, is very reassuring and gives us some confidence in our numerical scheme. 



3. Light-curves 

Once the jet dynamics and the pressure and density of the shocked fluid are known, 
the synchrotron plus inverse-Compton emissions are calculated from the fractional energies 
contained in the magnetic field and relativistic electrons which are parametrized by dimen- 
sionless numbers e B and e e respectively. The electrons are assumed to be accelerated to a 
power law distribution of energies, dN/d^ e oc % p , promptly behind the shock, and then cool 
due to radiative losses and adiabatic cooling. The local emissivity, j' u , (the energy emitted 
per unit volume per unit time, frequency & solid angle in the comoving frame), is approx- 
imated by a broken power law, with breaks at the cooling frequency, u' c , the synchrotron 
frequency, u' m , and the self absorption frequency, v' sa . 

Due to the curvature of the jet surface and its motion, photons arriving at some observed 
time t bs were emitted at different space-time points (r, t). The calculation of afterglow multi- 
wavelength light-curves takes into account appropriate integration over equal arrival time 
surface as given by the following equation for the flux density 

J?(t *W (1 + ^) 2 d 4 xf u ,(r,t) 

d L T 2 (r,t) [1 - v • n] 

where prime denotes quantities in the comoving frame of the fluid, n is the direction to the 
observer (in the lab frame), v' is related to the observed frequency v by the Doppler shift 
relation i.e. v' = v(\ + 2)7(1 — v • n), z and di are the redshift and luminosity distance 
of the burst, and d 4 x is the Lorentz invariant 4- volume element. The observer time t obs 
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is related to the lab-frame time t and location of the source r by: t b s = t — r • n, with 
r • h = r (cos 9 cos 9 ohs + sin 9 sin obs cos 0) . 

The lightcurves for the four jet profiles described in §2.3 are shown in Figure 5, and the 

micro-physics parameters for the shocked gas can be found in the figure caption. The average 

thermal Lorentz factor of shock heated protons is given by V th (9, t) = 1+711(0, t)/\u s {0, t) (7— 
2 

1)] — Il/(r fj, s ), which turns out to be in very good agreement with the value one obtains 
from shock jump conditions. 

The light curves for the Gaussian profile are rather similar to those for a 'top hat' jet, 
although the jet break for 'on axis' observers (0 o bs < 9 C ) is somewhat less sharp for the 
Gaussian profile, compared to a 'top hat' jet. The jet break is however still sufficiently sharp 
to be consistent with afterglow observations, in most cases. 

For a = and 6 = 2, the energy per unit solid angle (e) is independent of 9, so that after 
the deceleration time, tdec, the light curves for viewing angles, o bs < $dcc, become the same 
as for a spherical flow with energy Aire; note that for this jet model tdec T~ 8 / 3 oc 9^ s 3 and 
9 dcc oc t^ 6 . These light curves do not show a jet break, and are therefore not compatible 
with afterglow lightcurves that have a break. 

For a = 2 and 6 = (i.e. e oc @~ 2 and T =const, at initial time), there is a clear 
jet break for # b s > 29 c . The reason for this is that in this case much of the observed flux 
comes from small angles around the line of sight and the sharp break results when T<9~^ s . 
For # bs > 79 c there is a flattening of the light curve just before the jet break, due to the 
contribution from the inner parts of the jet, which is not seen in the observational data. 
This feature provides some constraint on this jet model and is discussed in Granot & Kumar 
(2002). 

For a = b = 2 (i.e. e, T oc Q~ 2 initially) the temporal decay slope before the jet break 
is steeper for small # bs and more moderate for large # bs, and the magnitude of the increase 
in the slope after the steepning of the lightcurve is larger for larger # bs- If such a correlation 
is found in the data, it would provide support for this jet profile. The jet break is smoother 
for small o b s , and sharper for large 9 ^ s , making it difficult to explain the sharp jet breaks 
and small inferred # Q bs (or opening angles for 'top hat' jets) that have been observed in quite 
a few afterglows. 

It can be seen that for a = b = 2 and a small r (~200 or less), the deceleration time, 
tdec, is quite large for large viewing angles, so that the light curve has a rising part at t < td ec 
(fig. 5). The fact that this is not seen in afterglow observations provides a lower limit on 
To or an upper limits on b, and thus can be used to constrain the structure of the jet. A 
more detailed analysis of the constraints that can be put on the jet profile from comparison 
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to afterglow observations is discussed in an accompanying paper (Granot & Kumar 2002). 

The light-curves obtained for the hydro simulation of jets can be reproduced, quantita- 
tively by one or both of two simple and extreme models, that are described below. For more 
details on these models and their results we refer the reader to Granot & Kumar (2002). 
The two different simple models are referred to as model 1 and model 2. 

In model 1, the energy per unit solid angle is assumed to retain its initial distribution, 
e{9,t) = e(9,t ). This represents the limiting case where there is very little lateral transport 
of energy. Model 2 attempts to make the opposite assumption, that is, it assumes the max- 
imal averaging of e over the angle 9 that is consistent with causality. The latter is achieved 
by averaging e over its initial distribution, over the range in 9 out to which a sound wave 
could have propagated from the initial time to- These two extreme assumptions are designed 
to bracket the expected range of possible behaviors for lateral energy transport. They are 
therefore expected to roughly cover the range of observed flux which a more rigirous treat- 
ment of the jet dynamics should give. In this sense they serve to quantify the uncertainties 
in the jet dynamics and light curves. For both models, the Lorentz factor is determined by 
energy conservation i.e. 



We have carried out hydrodynamical simulations of a relativistic, collimated, axisym- 
metric outflow propagating into an external medium. For simplicity, we used a uniform 
density medium for the calculations presented in this work. However, the numerical scheme 
developed in this paper is good for any axially symmetric external density distribution (in- 
cluding power laws with the distance from the source, etc.). We have reduced the problem 
to a 1-D system of partial differential equations by integrating over the radial thickness of 
the outflow, at a fixed lab-frame time, thereby greatly reducing the computation time. The 
hydrodynamical results were used to calculate the synchrotron emission and lightcurves for 
a variety of observer angles w.r.t. the symmetry axis. 

The model for GRB jets that is described in this paper is both rigorous and requires 
a very reasonable computational time, thus making it a useful tool for the study of GRB 
afterglows. In particular, it can be used in fits to afterglow observations, which can help 
constrain the jet profile, the external density profile, and the micro-physics parameters of 




(24) 



4. 



Conclusion 
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collisionless relativistic shocks. 

We find that for jets with smoothly varying energy per unit solid angle and Lorentz 
factor, T, the maximum transverse fluid velocity in the comoving frame of the shocked fluid 
is typically substantially less than the speed of sound: the peak velocity is of order 1/(T 89) , 
where 59 is the angular scale for the variation of the energy per unit solid angle, or of T. 
Thus large transverse velocities, approaching the sound speed, are realized only when the 
energy density varies rapidly with 9 or V decreases to (59)~ 1 . For the jet profiles examined in 
this paper, the largest lateral velocity occurred for the Gaussian profile, for which the initial 
gradients were largest. In fact, in this case a shock wave in the lateral direction developes 
due to the steep initial angular profile. This did not occur for the other jet profile which had 
a more moderate angular dependence. 

The energy per unit solid angle in the jet, e, is found to change very slowly with time 
for all four of the jet models we have analyzed, and to a reasonable approximation e can be 
considered essentially constant in time until T, along the jet axis, has dropped to ~ 4, which 
corresponds to about a week since the explosion in the observer frame (this is also the time 
when the transverse velocity is getting to be of order 0.2 c). 

Therefore, a simple model where the energy per unit solid angle is taken to be time 
independent, and each element of the jet behaves as if it were part of a spherical flow with 
the same e, can serve as a useful approximation for the jet dynamics, as long as the jet 
is sufficiently relativistic. As can be seen from Figure 6, this simple model indeed seems 
to reproduce the light-curves obtained using the hydro dynamical simulations quite well. 
However, the simple model must necessarily breakdown when the transverse velocity becomes 
of order c, and the energy density is no longer constant. Unfortunately, this is also the regime 
when our hydrodynamical calculation becomes unstable. 

We have calculated the observed light-curves in the R-band for several different jet 
angular profiles and viewing angles. We find that the light curves for a Gaussian jet profile 
are similar to those for a 'top hat' jet (as expected), and compatible with most observed jet 
breaks. The light curves for a jet with a constant energy per unit solid angle, e, but r(£ ) 
decreasing with 9, are similar to those for a spherical explosion, and thus not applicable to 
cases where we see a jet break in the light curve. For Jet profiles where initially e oc 0~ 2 
and T is either constant or oc 0~ 2 do produce jet breaks in the light curves, and have the 
advantage that these models can reproduce the 'observed' narrow range for the total energy 
in GRB relativistic outflows (Panaitescu & Kumar, 2002; Piran et al. 2001). A companion 
paper (Granot & Kumar, 2002) discusses some preliminary constraints on jet profiles from 
qualitative comparison with observations. 
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Fig. 1. — The different panels show the evolution of the Lorentz factor, T(0, t), the energy per 
unit solid angle, e(9, t), and the transverse velocity in the comoving frame, v' g /c = ug, for a jet 
with an initially Gaussian profile, i.e. T & e proportional to exp (-6 2 /26l) at the initial time. 
The parameters are: 9 C = 0.035 radian, e(0 = 0,t ) = 10 53 /4tt erg/sr, T{6 = 0,t ) = 200, 
and the density of the external medium is 10 particles per CC. At angles larger than about 
6.5° when T drops below 1.1, the LF and e are taken to be independent of 9. Note that the 
energy per unit solid angle does not change much with time except at large angles where it 
was small initially, and we see it increase with time. The sideways expansion of the jet can 
be seen in the bottom panel, which shows the edge of the "relativistic jet" , i.e. the vertical 
line where the transverse velocity drops to zero, which is moving to larger angles with time; 
the jet edge, and its motion, can also be seen in the top left panel as a sharp drop in Y — 1. 
The sharp jump in the lateral velocity and energy density may be understood as a formation 
of a shock wave in the lateral direction. The numerical scheme we use becomes unstable 
when T drops below 4 at 9 = or vg/c becomes larger than about 0.4; the LF near the 
edge of the jet is close to 1.01 when the code becomes unstable. It should be noted that the 
transverse velocity depends on the gradient of the energy density in the transverse direction 
and on the value of LF locally and not along the jet axis. And so our result for the transverse 
velocity - that it remaines below the sound speed throughout much of the jet evolution - is 
not compromised by the numerical instability. The lab-frame-time increases monotonically 
from the curve with smallest vg/c to the curve with highest vg/c (the lower panel), and the 
time increases from the highest to the lowest curves in the top two panels. 
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Fig. 2. — The different panels show the evolution of the Lorentz factor T(9,t), the energy 
per unit solid angle, e(9,t), and the transverse velocity in the comoving frame, v' /c = ug, 
for a jet with initial f = 1 + 199/(1 + 9 2 /9 2 c ) and e = e (l + 9 2 /9 2 c )- 1 with e = 10 53 /4vr 
erg/sr, and 9 C = 0.02 radian; this initial profile corresponds toa = 2& 6 = 2in the notation 
of equation (19). The density of the external medium is 10 particles per CC. Note that the 
memory of the initial core angle, 9 C , is not erased with time, and in fact it remains unchanged 
until quite late times (see top left panel). Moreover, the energy per unit solid angle does not 
change much with time either, except at large angles where it was initially small, and we 
see it increase slightly with time. The transverse velocity v' g (see lower panel) remains quite 
small throughout much of the evolution of jet and even at late stages is about a factor of 3 
smaller than the sound speed; the oscillations seen at early time in vg are almost certainly 
unphysical and numerical in origin. 
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Fig. 3. — Same as figure 2, expect that a = 2 & b = i.e. initially e = e (l + 9 2 /9 2 )~ 1 and 
T(9,t ) = 200, with 9 C = 0.02 radian. See caption for fig. 2 for other details. The minimum 
in vg at early times (the two lower curves on the bottom panel) is because v g is changing sign 
from negative values at small angles to positive values at larger angles and we have plotted 
\vg\. The reason for the negative velocity near the jet axis at early times is that in this model 
the deceleration time is largest at the pole (Rdec, ^dcc oc 0~ 2 / 3 = (1 + 9 2 /9 2 )^ 1 ^); thus the 
pressure integrated over the shell thickness is close-to- zero at the pole and non-zero at large 
9 at times much less that the deceleration timescale at the pole. In this case a pole-ward 
transverse flow ensues. At 9 > 9d cc (t) the e oc 0~ 2 profile dominates and induces a flow 
toward larger 9. At later time, much greater the deceleration time at the pole, the pressure 
integrated over the shell thickness, is indeed largest at the pole, as the referee mentioned, and 
the transverse velocity becomes positive everywhere. For (a,b) = (2,l) model the transverse 
velocity is everywhere positive at all times as expected from this argument and we have 
varified this numerically. 
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Fig. 4. — Same as Figure 3, expect that a = & b = 2, i.e. initially e(9, to) = 10 53 /47r erg/sr, 
and f = 1 + 199/(1 + 9 2 /9*), with 9 C = 0.02 rad. See caption of Fig. 2 for other details. 
Note that the jet becomes increasingly spherically symmetric with time in this case; T is 
independent of 9 up to an angle where shell deceleration has occurred, as in fact we expect 
when e is independent of angle. 
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Fig. 5. — Observed R-band light-curves for different viewing angles, obs , w.r.t. the jet axis 
are shown for the four jet profiles presented in Figs. 1-4; the basic jet model parameters 
are shown in the top right corner of each panel (see equation 19 for the definition of a 
& b). The remaining model parameters are the same for all panels: z — 1, the energy 
fraction in electrons and magnetic field are respectively e e = 0.5 & eg = 10~ 4 , the power-law 
index for electron distribution p = 2.5, and a constant external density n — 10 cm -3 . The 
lightcurves were calculated using the hydrodynamic simulation results shown in figures 1-4 
(see appropriate figure captions for the details of the jet models) and equation (23). For 
a = b = 2 the temporal decay slope before the jet break is more moderate for large 9 Q ^ S , and 
if T is not very large (200 in this figure) then the deceleration time for large 6> obs is rather 
large. For a = 2 and 6 = the jet break is rather sharp (as seen in observations), but at 
#obs ^ a few 9 C there is a flattening of the light curve just before the jet break (which is not 
seen in observations). 
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Fig. 6. — A comparison of lightcurves obtained with hydro simulation of jet (solid line) and 
a simple jet evolution model (dotted-curve) where the energy per unit solid angle is assumed 
to be time independent. The initial jet model is a power-law profile with a = 2, b = 2 and 
9 C = 1.1° (see equation 19 for the definition of a, b & 9 C ), and the observer location w.r.t. 
the jet axis, 9 f, s , is given in each panel. All the other parameters for the calculation are 
same as in figure 5 (see the caption for details). For clarity the observed flux for the simple 
model has been multiplied by a factor 2. 



